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LOW-MACH-NUMBER SLENDERNESS LIMIT 
FOR ELASTIC COSSERAT RODS 

FRANZISKA BAUS^, AXEL KLARi, NICOLE MARHEINEKE^’* *, AND RAIMUND WEGENER^ 


Abstract. This paper deals with the relation of the dynamic elastic Cosserat rod model and 
the Kirchhoff beam equations. We show that the Kirchhoff beam without angular inertia is the 
asymptotic limit of the Cosserat rod, as the slenderness parameter (ratio between rod diameter and 
length) and the Mach number (ratio between rod velocity and typical speed of sound) approach 
zero, i.e. low-Mach-number-slenderness limit. The asymptotic framework is exact up to fourth 
order in the small parameter and reveals a mathematical structure that allows a uniform handling 
of the transition regime between the models. To investigate this regime numerically, we apply a 
scheme that is based on a spatial Gauss-Legendre collocation and an Q-method in time. 


Keywords, dynamic elastic Cosserat rod, Kirchhoff beam, low-Mach-number-slenderness limit, 
asymptotic analysis, asymptotic-preserving scheme 
AMS-Classification. 74K10, 35Q74, 35B40, 35C20, 65Mxx 

1. Introduction 

The elastic rod theory is an old, extensively studied, but still current topic of research. Its 
foundations date back to, among others, Bernoulli, Kirchhoff [2^, the brothers Cosserat [12] and 
Love [28j . A comprehensive overview from today’s perspective is given in, for example [2j 136] . 
The investigations range from analytical aspects such as solution theory, stability and Hamiltonian 
structure to numerical methods (geometrically exact approaches, exact energy-momentum conserv¬ 
ing algorithms and symplectic schemes), of the broad literature see e.g. [29l [SO] HU HI [IHl HI] 
and [371 SOI HUHilll H]- Equally large is the field of applications: engineering mechanics (truss, 
fiber-reinforced materials, non-woven textiles, paper), biomolecular science (DNA, bacterial fibers), 
computer graphics etc. 

In this paper the specific focus lies on the asymptotic investigation of the relation between the 
dynamic Cosserat rod model and the Kirchhoff beam equations to describe the motion of slender 
elastic bodies. A rod in the special Cosserat theory [2] is represented by two constitutive elements: 
a parameterized time-dependent curve r and an attached orthonormal director triad {di,d2,d3} 
that specify the position and the orientation of the material cross-sections. Its deformation is due 
to tension, shear, bending and torsion. For general elastic material laws we derive beam equations 
of Kirchhoff-type in an asymptotic analysis, as the slenderness parameter e (ratio between rod 
diameter and length) and the Mach number Ma (ratio between rod velocity and typical speed of 
sound) approach zero. In the combined asymptotic limit e 0, Ma ^ 0 with Ma/e = n = const to 
which we refer as low-Mach-number-slenderness limit, the model equations show two characteristic 
changes: 1) the contact force becomes a variable to a constraint on the kinematics and the respective 
material law decouples from the model; 2) the angular inertia terms vanish. The terminology 
Kirchhoff stands for an inextensible and unshearable rod. Classically, the Kirchhoff constraint 
relates the beam tangent and the director triad, it is dsV = ds, see e.g. ladniEi]. As consequence 
of the constraint, the contact force and the strain change their roles in the model system. Their 
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proportionality operator is determined by the material law for the contact force, it can be identified 
with the moduli of linear elasticity theory. But in contrast to the linear theory whose validity is 
restricted to small deformations, large deformations are allowed as the material law can be chosen 
arbitrarily and the proportionality operator is just its linearization. Since namings of different beam 
models is frequently problematic and inconsistent in literature, we point out that we call the limit 
model a Kirchhoff beam when the system is equipped with the Euler-Bernoulli relation for the 
contact couple. In that case the deduced limit model that has additionally no angular inertia terms 
is also known as Kirchhoff-Love equations [23] . 

Deriving beam models with different degrees of freedom from the three-dimensional theory of 
elasticity is topic in several works, see for example the asymptotic limits in |101 m- The low- 
Mach-number-slenderness asymptotics presented in this paper reveals a direct scaling between the 
Cosserat rod and the Kirchhoff beam without angular inertia and complements the previous works. 
The result has its analogue in the fluid dynamics where the low-Mach-number asymptotics describes 
the transition from a compressible to an incompressible fluid EZlilllSS]- For elastic rods the speed 
of sound is not a unique quantity, instead compression and shear waves induce different speeds of 
sound that imply different Mach numbers. For the asymptotic derivation, it is assumed that there 
exists a typical magnitude of the contact force. Such a scaling presupposes that the Mach numbers 
associated with compression and shear behave similar. This corresponds particularly to a constant 
Poisson ratio for a linear isotropic material, as it is often considered in macroscopic applications. 

The asymptotic rod behavior is here characterized by the slenderness parameter e and the (typi¬ 
cal) Mach number Ma which we relate according to Ma/e = /i = const. In the transition regime for 
small e (e —>• 0) the type of the model equations changes: time derivatives degenerate, the system 
becomes stiff with a constraint. The asymptotic framework reveals a mathematical structure that 
allows a uniform robust numerical handling of this regime. We show that the asymptotic systems 
(limit system and its first-order correction) which follow from an even power series expansion are 
together exact up to order 0{e*). Moreover, the conservation of energy is ensured, if hyper-elastic 
material laws, external potential forces and appropriate boundary conditions are presupposed. Hav¬ 
ing a Newton method for the numerical treatment in mind, solving both asymptotic systems is of 
similar computational effort as solving the original e-dependent system while it is much more ac¬ 
curate and robust in computing the influence of small e-values. We point out that the underlying 
discretization is replaceable. In this paper we propose an asymptotic-preserving scheme in the spirit 
of the generalized a-method that can be switched between energy-conserving and dissipative [gEg. 
This is advantageous for applications with external non-potential forces, such as fiber-fluid interac¬ 
tions in non-woven manufacturing [saEi] or hair simulations |4| [42] . We refer to |39l [35] for exact 
energy-momentum conserving algorithms and to HIS] for discrete geometric approaches. More¬ 
over, we refer to the research work in fluid dynamics where the low-Mach-number asymptotics has 
been used to develop and extend numerical schemes for the compressible-incompressible transition 
regime, e.g. mill]. 

This paper is structured as follows. Starting with a brief introduction into the special Cosserat 
rod theory, we present the low-Mach-number-slenderness asymptotics in Section |2] We derive 
the beam equations of Kirchhoff-type for general elastic material laws. To study numerically the 
performance of the asymptotic framework in the transition regime for small e, we consider a two- 
dimensional Fuler-Bernoulli cantilever beam in Section [3] Three variants of the underlying rod 
model can be distinguished, depending on the chosen kinematic/geometric formulation. They imply 
temporal or spatial differential-algebraic systems of different index after semi-discretization in space 
or time, respectively. In the Appendix we provide details to the underlying numerical scheme that is 
applicable to all formulations in the transition regime. It is based on a Gauss-Legendre collocation 
in space (finite differences) and an a-method in time, but can be also viewed as a conservative 
finite-volume method. 
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2. Asymptotic relation between elastic Cosserat rod and Kirchhoff beam 

2.1. Special Cosserat rod theory. An elastic thread is a slender long body, i.e. a rod in three- 
dimensional continuum mechanics. Because of its geometry the dynamics might be reduced to a 
one-dimensional description by averaging the underlying balance laws over its cross-sections. This 
procedure is based on the assumption that the displacement field in each cross-section can be 
expressed in terms of a finite number of vector- and tensor-valued quantities. The special Cosserat 
rod theory [5] consists of only two constitutive elements in the three-dimensional Euclidean space , 
a curve r : —>■ specifying the position and an orthonormal director triad {di, d 2 , da} : V 

characterizing the orientation of the cross-sections. In S | s S [sojSb], t > 0}, the 

parameter s denotes the material cross-section (material point) and t the time. The rod system 
involves four kinematic/geometric and two dynamic equations (balance laws for linear and angular 
momentum) 

dtr = V, dsV = T 

dtdk = w X dk, a^dk =/« X dk, fc = l,2,3 

dt{{pA)v)=dsn + i 
5t((pJ) • u}) = dsHi -I- X X n -I- 1 

with tangent r, generalized curvature k, linear velocity v and angular velocity U3. The mass line 
density (pA) is time-independent for materially closed systems, but modeled as time-dependent 
for applications, such as evaporation and aggregation. The tensor-valued moment of inertia of the 
cross-sections (pJ) depends on the configuration of the triad and is hence always time-dependent. 
To close the system of equations we need to specify the external loads f and 1, boundary and initial 
conditions as well as elastic material laws for the contact force n and couple m. In this paper we 
consider time-independent mass properties of the rod, i.e. dt{pA) = 0 and 9t(di • (pj) ■ dj) = 0 for 
i,j = I, 2, 3, and neglect external couples, i.e. 1 = 0, for reasons of a simple presentation. However, 
extensions are straightforward possible. Alternatively to (EtD, the rod system can be also set up 
by using the compatibility relations between the kinematics and geometry for the curve and for the 
triad. These two compatibility conditions 

dtT = dsV 

dtn = dgU) + u X K 

replace then either the two kinematic equations or the two geometric equations. 

Notation 1 (Model variants (M), (T), (S)). Depending on the chosen kinematic/geometric formu¬ 
lation, we distinguish three variants of the rod model: 

(M) kinematic and geometric eguations and balance laws (cf. (12.11 ) ) 

(T) kinematic equations, compatibility conditions and balance laws 
(S) geometric equations, compatibility conditions and balance laws 
Throughout the paper we often summarize the three systems in a compact form (M-T-S) for read¬ 
ability. Apart from the balance laws, (M-T-S) contains then all six kinematic, geometric and com¬ 
patibility equations together, see for the first time in (1221). 

The model variant (T) yields a hyperbolic system for a dynamic elastic rod, whereas (S) is very 
suitable for the transition to a stationary consideration. Presupposing hyper-elastic constitutive 
relations and external potential loads, the variant (M) is known as bi-Hamiltonian form in literature 
[Si US]. The choice of the formulation can affect the numerical simulations as we will comment on 
in Section li and the Appendix. However, the distinction plays no role for the asymptotics, thus we 
make use of the compact notation (M-T-S) in the following. 

For the objective formulation of the material laws it is useful to rewrite the rod model in the 
director basis. To an arbitrary vector field z = X]i=i ^ indicate the 

coordinate tuples corresponding to the director basis {di,d 2 ,d 3 } and to a fixed outer Cartesian 
basis { 61 , 62 , 63 } by z = (zi, 22 , 2 : 3 ) S R^ and z = ( 21 , 22 ,- 23 ) S R^, respectively. The director 
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basis can be transformed into the outer basis by the tensor-valued rotation D, i.e. D = ej (g) di = 
DijBi g) ej £ g with the associated orthogonal matrix D = (Dy) = (dj • ej) £ iS'0(3). For the 
coordinates, the relation D-z = z holds - as well as D-c^tz = dtz + oj x z and D-dgZ = OsZ + k x z. The 
rotation matrix D can be parametrized, for example, in Euler angles or unit quaternions. Moreover, 
canonical basis vectors in are denoted by e^, i = 1,2, 3, e.g. ei = (1, 0,0). In the director basis, 
the Cosserat rod model has the following form (M-T-S) 

D • dtf = V, D • dgf = r 

dtD = —cli X D, dsD = —K X D 

dtT = ds'J + KX\/—UJXT 

dtK = dsUj + K X uj ( 2 . 2 ) 


{pA)dt\/ = dsn + K X n — UJ X (pv4)v + D • f 
(pj) ■ dtuj = cism -l-KXm-l-rxn—wx ((pJ) • uj) 
with general elastic material laws 

n(s,t) = J\f{T{s,t),K{s,t),s), m{s,t) = M{T{s,t), K{s,t), s). 

Following the assumption on the mass properties, (pJ) is here the time-independent representation 
of the inertia tensor in the director basis. In the stated general form the objective elastic constitutive 
functions Af and A4 depend on the strain variables (coordinate tuples) r and k whose components 
measure shear ri, T 2 , dilatation T 3 , flexure ki, K 2 and torsion K 3 . The subsequent asymptotic 
analysis is valid for material laws that possess the following properties. 

Assumption 2 (Properties of the material law A/”). Assume that the material law Af for the contact 
forces fulfills the following conditions: 

a) Af{T, K,s) = 0 holds if and only if t = r(s). 

b) Af is sufficiently regular and drAf{T{s),K,s) is an invertible linear operator (matrix) for 
arbitrary k and s. 

Remark 3. Assumption\^) presumes that the family of stress-free configurations is uniquely deter¬ 
mined by a strain field f. Without loss of generality an arc-length parameterization can be chosen for 
these configurations, implying HfH = 1. Assumption\^) is naturally satisfied by the monotonicity 
condition on the constitutive laws [2]. In case of hyper-elastic constitutive relations, i.e. Af = 
and A4 = d^'I', it is expressed in terms of a strictly convex elastic potential 'b. 


Example 4 (Specification of material laws). As an example for a (hyper-)elastic constitutive law 
satisfying the assumptions above one can think of a Timoshenko beam under small deformations 
that is equipped with the affine linear Euler-Bernoulli relation for the contact couple, i.e. 

Af{T, K, s) = (EA)(s) • (r - f(s)), A4{t, k, s) = (EJ)(s) • (k - k{s)). 


In literature the Timoshenko beam is in general associated with r(s) = 63. The positive-definite 
tensor-valued functions (EA) and (EJ) can be expressed by the scalar-valued properties of the cross- 
section associated Young’s modulus (EA) and the shearing modulus {GA) - in combination with 
the mass line density {pA) and the moment of inertia {pi). In case of a geometry with circular 
cross-sections, we obtain the diagonal /orm|3 


(pJ) = (p/)P2, 


(EA) 


(EA) ^ 

a 


(EJ) 


{PI){EA) 

ipA) 


P 2/a 


where Pfe = diag(I, l,k), fc £ R., and a = {EA)/{GA) = 2(1 -|- u) with Poisson ratio v £ [0,0.5). 
The physical quantities {pA), {pi) and {EA) are particularly constant for an homogeneous thread. 


*The diagonal forms of (pJ) and of the linear operator associated with A4 are necessary for the embedding into a 
2d test scenario. But note that for this purpose the special choice of circular cross-sections is not compulsory. 
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2.2. Asymptotic low-Mach-number—slenderness limit. Proceeding from the Cosserat rod 
model, we derive beam equations of Kirchhoff-type (a generalized string model) as the combined 
slenderness and low-Mach-number limits in this subsection. We consider an elastic thread. It has 
various physical properties for which we choose typical constant reference values that we mark by 
the subscript *, i.e. mass line density (pA)*, moment of inertia (pi)*, length L*, magnitude of 
mean velocity 14 and contact force A*. When making (12.2|) dimensionless, we can characterize the 
dynamics of the thread by help of two dimensionless numbers: the slenderness parameter e that is 
the ratio between the thread’s diameter and its length as well as the (typical) Mach number Ma 
that is the ratio between the thread’s velocity and the speed of sound 


1 


' {pi)* 

{pA)* 


Ma = W, 


'{pA). 


A* 


The speed of sound is here determined by the typical contact force A*. In the classical sense of 
compression and shear waves, an elastic thread has certainly different speeds of sounds. Hence, the 
applied scaling presupposes that the associated Mach numbers behave similar. This means in the 
special case of a linear isotropic material (cf. Example |4]) that the typical Young’s modulus could 
be chosen. A* = (AA)*, and that the Poisson ratio satisfies i/ = const in the asymptotics. 

For the scaling, we use the following reference values 

so = ro = L*, Vo = 14, to = Wo = kq = (pA)o = (pA)*. (pJ)o = {pi)* 

V* Li, Li, 


no = {pA)*V^, mo = {pA)*L*V^, No = A*, Mo = 


{pI)*N* 


fo = 


{pA)*V^ 


{pA)*L* 

and introduce to every dimensional variable z the associated dimensionless one z as 

z(sos, tot) = zoi{s, t), s = sqs, t = tot. 

Skipping the superscript “ for readability, the dimensionless rod system is then given by 
D • 9tf = V, D • i9sf = T 

dtD = — w X D, i9sD = —K X D 

dtT = ds'J + + 

dtK = dgOJ + K X oj (2.3) 


[pA)dtsi = dsn -I- K X n -|- (pA)v x w -|- D ■ f 
e^(pJ) • dtOJ = dsvn -|-KXm-|-rxn-|- e^((pJ) • w) x 
with general elastic material laws {N satisfying Assumption [5]) 


'^{^A) = T^Af{'r{s,t),K{s,t),s), 
Ma 


’^{S,t) = Z^M{Tis,t),K{s,t),s). 
Ma 


To establish an asymptotic relation between the Cosserat rod (loi) and a beam of Kirchhoff-type, 
we consider the limits 


0, Ma —t 0 


with 


Ma 


= p, p = const > 0 


(i.e. asymptotic limit along straight lines in the (e, Ma)-plane with slope p). As the small parameter 
e appears only in second power in the equations (e^, Ma"^ = p^e^), we expand all quantities in an 
even power series with respect to the slenderness parameter, i.e. z = e^*^*^*^ The contact force 

n{s,t) = {p^)~^Af{T{s,t),K{s,t),s) implies 


0 = N{t^°\s, t), t), s), as e = 0, 


c(o) 


hence r^°^(s,t) = t{s) 
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according to Assumption • Consequently, we set the strain to be r = r + with a function 
^ ~ 0(1). Then the contact force becomes n(s,t) = {t{s) + e^^(s, t), k(s, <), s), whose 

leading order we obtain via Taylor expansion 

{drN'{T(s),K^°\s,t),s) + d^Af iT(s), K^°\s,t), s) • 

Here, the second term vanishes since M(t(s), k, s) = 0 for all k and s (Assumption [2]). So, the 
expression ^ ■ n with L = drAf{T(-), k, ■) is exact up to an error of O(e^). Inserting this 

expression for ^ into (12.31) yields consequently a e-dependent consistent system up to order O(e^) 

D • 9tf = V, D • dgl = r -I- ■ n 

dtD = —Lo X D, 5sD = —K X D 

• dtn = dgV + KXy + TXU} + • n x w 

dtK = dsOJ + KX UJ (^■^) 


with 


{pA)dty = dsn + K X n-f (pA)y x w -I- D • f 
e^(pJ) ■ dtuj = 9sm -|-KXm-|-'fxn-|- ■ n x n -|- (pj) ■ oj x oj) 


vn(s,t) = P ^A4(t(s) + e^p^L ^ • n(s, t), k(s, t), s), L(s, t) = 5 t-A/’(t-(s), n(s, t), s). 

Setting e = 0 in (12.4|) . we obtain the low-Mach-number-slenderness limit for elastic bodies. The limit 
equations (12.Sp describe a simplified model of the special Cosserat rod theory, where the angular 
inertia terms vanish and instead of a material law for the contact force the time-independence of 
the strain held is imposed (dtT = 0). This property enforces the limit beam to be inextensible and 
unshearable and is known as Kirchhoff constraint. Classically, the Kirchhoff constraint relates the 
tangent and the director triad r = da, it is often stated as D • dgf = r = 63 in the director basis, see 
e.g. unniiis]. This corresponds to a straight rod curve with perpendicular cross-sections as choice 
for a stress-free conhguration, f = 63 . 

D • 5tf = V, D • 5sf = T 

dtD = —OJ X D, 5sD = —K X D 

0 = 5sV -l-Atxv-l-rxtLi 

dtK = dgUJ + K X OJ (^■^) 


(pA)dty = 5sn -|- k x n -|- (pA)y x w -I- D • f 
0 = 5sm-|-KXm-ffxn 


with 

m(s,t) = p~'^M(t(s), K(s,t), s). 

In the formulation of (12.41) the contact force n and the strain r change their roles, their proportion¬ 
ality operator L is thereby determined by the material law Af. In the limit e = 0 in (1^ . n becomes 
a variable to the constraint 5(T = 0, whereas the material law decouples and is not longer needed 
for the determination of the solution. The structure of (12.41) obviously changes from a hyperbolic to 
a degenerate differential-hyperbolic-like system with constraint as e —>■ 0, the system becomes stiff. 
Note that L can here be identihed with the moduli of the linear elasticity theory, linear operator 
(EA), cf. Example m But in contrast to the Timoshenko beam whose validity is restricted to small 
deformations, system (12.41) allows for large deformation as Af can be chosen arbitrarily and L is just 
its linearization, L = drAf{T{-), ■). Without loss of generality, we might thus restrict our consider¬ 
ations for small e to an (affine) linear material law of the form Af{T, k, s) = L(s) • (r — 'f(s)) which 
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implies an explicit and easily invertible relation between contact force and strain. This relation 
offers advantages in setting up a numerical scheme that is applicable to both, limit and e-dependent 
system as e —>■ 0 . 

The presented low-Mach-number-slenderness limit finds its analogue in fluid dynamics where 
the low-Mach-number limit represents the incompressible limit for compressible fluids Here, 

in elasticity, the asymptotic derivation goes with slow dynamics and slenderness. The limit is 
valid for arbitrary elastic material laws A4 for the contact couple. When the limit system (12.5p 
is equipped with the Euler-Bernoulli relation it describes a Kirchhoff beam. In that case with 
f = 63 the system is also known as Kirchhoff-Love equations [23]. However, in the classical theory 
the terminology Kirchhoff beam is much wider and stands for an inextensible and unshearable rod 
(kinetic analogue |2]). The vanishing of the angular inertia terms is generally not presupposed for 
a Kirchhoff beam (see e.g. mM), but results here as consequence of the chosen scaling and the 
associated asymptotics. 

2.3. Asymptotic framework with Euler-Bernoulli material law. In the further work we 
focus on the asymptotic framework in a special case, where we restrict on a linear Euler-Bernoulli 
relation = (EJ) • k and on an homogeneous thread with circular cross-sections and f = 63 . Note 
that these are just technical simplifications to facilitate the numerical studies. They can be easily 
dropped if it is relevant for certain applications. 

In this case the dimensionless quantities become 

{pA) = l, (pJ) = P2, L = (EA)=a-iPa, (EJ) = P2/a. 

where P^ = diag(1,1,/c), k € M. and a = 2(1 -|- v) with Poisson ratio v € [0,0.5) - in accordance 
to Example |3| The physical quantities stated in Example Sj correspond to the reference ^-values 
chosen for the scaling. Hence, we deal with the following system in the (M-T-S) notation that need 
to be supplemented with appropriate, consistent initial and boundary conditions, it describes the 
Cosserat rod for e > 0 and the Kirchhoff beam for e = 0, cf. (ED), (USD 

D • 5tf = V, D • dsf = 63 - 1 - e^p^aPi/a ■ n 

5(D = —uj X D, 9s D = —K X D 

e^p^aPi/a • 9(11 = 9sV -|-KXv-|-e 3 Xa;-|- e^p'^a{Pi/a • n) x w 

dfK = dsUj + KXUJ 0 ) 

dt\/ = dgn + K X n-l-vxw-l-D-f 
e^P2 • 9(0; = fj,~^ (Pa/a -dsK + nx (P2/Q • k)) -E 63 x n 
-E (/i^a(Pi/a ■ n) X n -E (P2 • o;) X to) 

Remark 5 (Re-formulations of the limit system). In the limit system, e = 0, the Kirchhoff con¬ 
straint poses a geometric relation between curve and director triad in favor of a material law for 
the contact force. This motivates the so-called centerline-angle representation of the elastic Kirch¬ 
hoff beam [25j that renounces the evaluation of the director triad. Alternatively, the Kirchhoff-Love 
equations might be known in the invariant form of a wavelike equation for r with small elliptic 
regularization due to the bending stiffness (incorporated in p) jU] l3T] 

dttr = ds{Tdsr) - p~‘^dssssY-\-M{dsr X dsssr)-\-{, dsM = 0, ||9sr|| = 1 (2.7) 

with torsion couple M = K 3 . The modified traction T = n-da — ^“^||9ssr|P acts thereby as 

Lagrange multiplier to the constraint. Obviously, both re-formulations of the limit system reduce the 
number of variables, but strongly change the equations ’ structure. Appropriate numerical schemes 
can be found in the respective literature. However, through the low-Mach-number-slenderness limit 
that clarifies the asymptotic relation between (1^ and (EH) a uniform numerical treatment for 
e —>■ 0 JS made possible. 
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The asymptotic framework provides a special structure of the equations that is exploited for 
the numerical handling. Let 4> denote the vector-valued function comprising all system variables, 
0(s,<) € M’”, then the e-dependent partial differential algebraic system (12.61) can be summarized 
for all model variants (M), (T), (S) in the form 


A,-at0 + B,-a,cD + c,(cD) = 0, (2.8) 

where Aj and are - possibly singular - matrices with constant coefficients and Cg is a vector-valued 
nonlinear function in 0. Their e-dependence can be expressed as A^ = A^^^ -|- e^A*^^\ analogously 
for Bj and c^. Expanding 0 in the even power series 0 = X]r=o and plugging it in (12.81) - as 

before we find the Kirchhoff beam in 0(1) and its first-order correction in O(e^) 

A(°) ■5t0(°)-I-B(°) •5s0(°)-f c(°)(0(°)) =0, (2.9) 

A(°) • dM^'> + B(°) • 9,0(1) -f 90c(°)(0(°)) ■ 0(1) = f[0(°)], (2.10) 

with - f[0°] = A(i) • 9(0(1) -f B(i) • 9,0(1) -f c(i)(0(i)) 

where 9<tc(i) denotes the Jacobi matrix of c(i). Both asymptotic systems together (I2.9I) - (I2.10I) 
are exact up to order 0{e'^). As usual for asymptotic considerations, the systems have a similar 
equation structure with the same system matrices A(i) and B(i), they just differ in the right-hand- 
side and the dependence on the variable. The first-order correction (12.101) is trivially linear in the 
variable, whereas the limit system (12.91) keeps the nonlinearity of the original e-dependent system 
(ITSI) . Having a Newton method for the numerical treatment of the nonlinear term c(i)(0(i)) in 
mind, we need its Jacobi matrix already for the computation of the limit system. So, the assembly 
of the linear system matrix associated with the first-order correction is for free. 

Remark 6 (Conservation of energy). The Cosserat rod model as well as the Kirchhoff beam equa¬ 
tions in (1^ are energy-conserving for external potential forces |38j . presupposing classical bound¬ 
ary conditions such as, for example, a cantilever beam or a beam with stress-free ends. This holds 
also true for the system associated with the first-order correction. The corresponding energy with 
(external) potential energy V is given by w, = -f e^w(i) -f O(e^) 

u;(i) = i(v(0))2 + ^(.(1) . Py, . (,(1) + E(f(i)) 

((;(!) = v(i) ■ v(i) + ^n(i) • P 2 /a • «(1) + ia;(i) • P 2 • c.(i) + ^n(i) • P^/, • nd) + ^(fd)). 


3. Numerical studies 

Being interested in the numerical performance of the asymptotic framework as e —> 0, we study 
and discuss the asymptotic convergence and efficiency of the approach in this section. To make use 
of the structure of the asymptotic systems we apply a Newton method to the discretized equations. 
The discretization is replaceable. We apply here an asymptotic-preserving scheme in the spirit of the 
generalized a-method [giig to the underlying systems of partial differential algebraic equations. 
This scheme leaves the freedom to be switched between energy-conserving and dissipative. This 
feature can be advantageous when dealing with applications with external non-potential forces 
(such as fiber-fluid interactions in fiber spinning |31) , non-woven manufacturing EH, paper making 
[18j or hair simulations 0)- For details on the numerical method we refer to the Appendix. As test 
case we use the two-dimensional Euler-Bernoulli cantilever beam under gravity m- This test case 
offers the advantage that the number of model variables can be reduced from to = 19 to 9, while 
the e-dependent structure of the equations that is of interest is kept. So, clarity is given for the 
investigation. 
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3.1. Test case. Let { 01 , 62 , 63 } be the outer Cartesian basis. Consider a thread fixed at one end 
(s = 0 ) and stress-free at the other end (s = 1) that is initially static, stress-free and straight 
in direction of 61 . Let it be exposed to transversal oscillations in the 6 i-e 2 -plane due to gravity 
f = —{pA)ge-i such that d 2 = 63 during its motion (cf. Figure 13.11) . For this two-dimensional 
scenario, the model system (|2.6I) in the (M-T-S) notation simplifies to 

D(a) ■ dA = V, D(q!) • ds'r = 62 -I- e^/r^aPi/o • n 

dta = uj, dsa = k 

■ dtr\ = ds'J — — wei -|- e^/i^wPa • 

dtK = dsLO 


i9tv = Ssti — -|- D(a) • f 

e^dtOJ = p.~^dsK -I- ni -I- e^^^(l — a)nin3 


with the dimensionless force f = —Fr ^ 


62 and the respective initial and boundary conditions 


f(s, 0 ) = sei, q;(s, 0 ) = k(s, 0 ) = w(s, 0 ) = 0 , n(s, 0 ) = v(s, 0 ) = 0 

f( 0 ,t) = v( 0 ,t) = 0 , a( 0 ,t) = a;( 0 ,t) = 0 , n(l,t) = 0 , K(l,t) = 0 . 

Here, we use the abbreviations 


f=(ri,r 2 ), f = (/i,/ 2 ), v=(ui,U 3 ), n = ( 711 , 713 ), w = W 2 , k = K 2 , m = m 2 , 

Pfc = diag(l, k), k = R. Moreover, = (—Z 3 , zi) denotes the tuple perpendicular to z for z G {v, n}. 
The rotation matrix D is expressed in terms of the single angle a 


D(a) 


— sin a cos a 
cos a sin a 


a = Z( 6 i,d 3 ). 


The Froude number Fr represents the ratio of inertia and gravity. 

Because the asymptotic results turn out not to depend very sensitively on the specific choice of 
parameters, we use the setting (Fr,/i,a) = (1,10,2.5) with the end time T = 2.5 as example in 
the following studies. This parameter setting is characteristic in the context of non-woven manu¬ 
facturing [321 El]. Typical properties of polymer fibers are d = 10 ^m (diameter), L = 10 ^ m, 



Figure 3.1. Test case: dynamics of the 2d Euler-Bernoulli cantilever beam under 
gravity f. Simulation of the limit system (e = 0) and the e-dependent system 
(e = 0 . 02 ). 











10 


F. BAUS, A. KLAR, N. MARHEINEKE, AND R. WEGENER 


p = lO^kg/m^, E = 10^°kg/(ms^), u = 0.25 and V = Im/s, implying e ^ The end 

time T is chosen respectively with regard to stability issues (see Remark [7]). The dynamics of the 
cantilever beam under gravity is illustrated for e = 0 and e = 0.02 in Figure [5TT] 

Remark 7. For the planar inextensible beam an existence proof is derived in the absence of gravity 
and for non-vanishing angular inertia in [7]. Nonlinear planar and non-planar responses of the in¬ 
extensible beam were investigated numerically in [31] using a combination of the Galerkin procedure 
and the method of multiple scales, it turned out that the nonlinear geometric terms produce a hard¬ 
ening effect and dominate the non-planar responses for all modes. The inertia terms particularly 
dominate the high-freguency modes. A global bifurcation analysis for a cantilever beam subjected to 
a harmonic axial excitation and transverse excitations at the free end was performed in [44j . Taking 
into account these results, the simulations of our limit system that has no angular inertia are run 
in the stable planar regime. 

3.2. Results. In the asymptotic framework we have 4>e = 0^°^ + + O(e^), i.e. the limit 

system and its first-order correction are analytically exact up to the order cf. (I2.8I) - (I2.10I) . 

The solutions of the limit system 0^°^ (Kirchhoff beam), its first-order correction 0^^^ as well as of 
the original e-dependent system 0^ (Cosserat rod) are numerically approximated by and 

<p^. We observe that the numerical approximations inherit the asymptotic relation as desired, this 
means that 

= cl ~ 0(e2), cl =£2 Cl (3.1) 

(p^ - = C20(e‘^), C2=e‘^C2. (3.2) 

hold true. Figure 13^ shows the £^(0, l)-norm of the vector-valued functions c? and Ci, i = 1, 2, in 
dependence on e for e G [10“^°,1]. The quantities are computed at time T = 2 for the different 
model variants (M), (T) and (S) introduced in Notation |T] by applying the numerical scheme with 
the parameters At = As = 10“^ and A = 1. Comparing the model variants, they all yield the 
quadratic convergence for the first term cl with the same (e-independent) magnitude ||ci|| in the 
range e G its relative deviation from lies below 10“^. For the second term 

C 2 the model variants show the quartic convergence in the range e G [10“"^, 10“^] and (S) even in 
e G [10“®,10“^]. The reduced order of convergence (down to p = 2 in e G [10“®, 10“'*]) for (T) 
and (M) might be explained by the abandonment of boundary conditions. All boundary conditions 
are explicitly prescribed and incorporated in the scheme for (S), whereas the boundary conditions 
associated with f, a for (T) and to v, n for (M) follow only numerically from the initial conditions 
and the model equations. The approximations are consistent but the effect of the term in 

the numerical differentiation is below the computational accuracy (cancellation error). Also the 
tails that are observed for very small e in the convergence plots come from numerical noise. The 
striking switch in the convergence behavior of both terms c?, i = 1 , 2, that occurs at e = 10 “^ for all 
model variants can be explained by the asymptotics itself because the asymptotic framework holds 
true only for e sufficiently small. This can be also clearly seen in Figure [33] where the beam curve 
associated with the comparatively large e lies away from the limit curve. Note that the curves would 
be indistinguishable for e < 10 “^. The corresponding peaks observed at e = 10 “^ in Figure [331 
(right) leave freedom for speculations: they might be due to the transversal stress component whose 
effect is then superposed by the other variables. 

As the Cosserat rod model contains the small asymptotic parameter explicitly, the system of 
equations becomes stiff and difficult to solve numerically as e —>■ 0. The asymptotic framework 
provides a special structure of the asymptotic systems (I2.9I) - (I2.10I) that can be exploited in the 
numerics when solving the discretized equations with a Newton method. The nonlinear original e- 
dependent system as well as the nonlinear limit system require 2-3 Newton iterations in average per 
time step for all model variants, supposing At ^ 0(10“^). The resulting linear systems are of same 
size and have a similar block-structured band matrix. The Jacobian of the limit system equals the 
linear system matrix associated with the first-order correction. Therefore, the effort of computing 
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Figure 3.2. Conservation of the asymptotic relations for the different model vari¬ 
ants (cf. (13.Ill - (13.21) . = e^*Ci). Top: first term ||ci|| (left) with magnitude ||ci|| 

(right) plotted over e. The solid line with p = 2 indicates quadratic convergence. 
Bottom: second term Ijc^H (left) with magnitude ||c 2 || (right) plotted over e. The 
solid lines with p = A and p = 2 indicate quartic and quadratic convergence, re¬ 
spectively. 


the two asymptotic systems is similar to the effort for the original e-dependent system, as the first- 
order correction costs only a single linear solving. Moreover, the solving of the asymptotic systems 
is much more accurate and robust in determining the influence of small e-values. Consequently, 
the asymptotic framework makes a uniform numerical handling of the transition regime for small e 
easily possible. The choice of the model variant has only a very marginal effect on the numerical 
performance. It influences neither efficiency nor convergence properties, for details see the Appendix. 
However, (T) seems to be preferable for energy-conservation and (S) for robustness and asymptotic- 
conservation. 


4. Conclusion 

The low-Mach-number-slenderness limit (e —>■ 0 and Ma 0 with Ma/e = p. = const > 0) de¬ 
scribes the asymptotic relation between the dynamic elastic Cosserat rod and the inextensible, un- 
shearable Kirchhoff beam without angular inertia. Its derivation requires the unique determination 
of the family of stress-free configurations by a strain field and the monotonicity of the constitutive 
law for the contact force. The limit is valid for general material laws for the contact couple. As 
the Cosserat rod model contains the small asymptotic parameter explicitly, the system of equations 
becomes stiff and difficult to solve numerically as e —>■ 0. The derived asymptotic framework allows 
a uniform numerical treatment of the transition regime between the two models. We point out that 
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solving the asymptotic systems (limit system and its first-order correction) that are exact up to 
order 0 {e'^) is of same computational costs as solving the original e-dependent system, but much 
more robust since the equations are independent of e. Hence, we propose to use the asymptotic 
systems for the numerical investigation of applications where the beam velocity compared to the 
typical speed of sound is small and the slenderness ratio is small, as it is the case, for example, in 
fiber dynamic simulations of non-woven production processes EH 131] or hair modeling in computer 
graphics [Slij- 


Appendix A. Numerical Method 


The rod models consist of systems of partial differential algebraic equations. The model variants 
(M), (T) and (S) introduced in Notation [T] imply temporal or spatial differential-algebraic systems 
of different index (index 0 to 3) after semi-discretization in space or time, respectively. We propose 
a numerical scheme that is applicable to all formulations in the asymptotic framework. For the 
numerical treatment of the respective model equations we combine a Gauss-Legendre collocation 
method (finite differences) on a equidistant space grid with a flexible time integration. It can be 
also viewed as a conservative finite-volume method. Using a temporal tuning parameter A G [0.5,1] 
the time integration can be switched continuously from an energy-conserving Gauss midpoint rule 
(A = 0.5) to a dissipative implicit Euler method (A = 1) in the spirit of a generalized a-method 
The discretized Cosserat rod and Kirchhoff beam models result in nonlinear systems of 
equations that are solved with a Newton method with Armijo step size control. 


Remark 8. To exploit the structure of the asymptotic systems only the application of the Newton 
method is essential, whereas the underlying discretization is replaceable. Hence, we refer to the 
well-established results by Simo & Vu-Quoc, Simo et al. gnisiiiss] or Romero & Armero [35] for 
exact energy-momentum conserving algorithms. Preserving conservation properties and objectivity 
is also the topic in, among others, where the beam equations are formulated as a constraint 

Hamilitonian system. The ingredients of the scheme are finite elements and a G-equivariant discrete 
derivative, the constraints are treated by the Lagrange multiplier method, the penalty method or the 
augmented Lagrange method. Ln [24] the Cosserat rod model is transformed into a Lagrangian 
differential-algebraic system of index 3 that is reduced to index 0 by introducing Baumgarte penalty 
accelerations for stabilization and solved by using finite differences on a staggered grid. The stiffness 
of the inextensible, unshearable beam with angular inertia is handled by the schemes in e.g. [HE]. 
Further approaches based on spatial semi-discretizations and Lagrangian mechanics can be found in 
literature; for algorithms used in the computer graphics see e.g. [10I33]- 


A.l. Numerical scheme. Consider the spatial and temporal grids Si = iAs, i = 0,...,N and 
P = jAt, j = 0,1,... with hxed cell and step size As and At. Let the vector of the unknown system 
variables at the grid point be denoted by ifl G K™ and all spatial information at the time 

level P be summarized in ip^ G We treat (si+i/ 2 ,P~^^) with A G [0.5,1] as collocation 

points. The idea of the scheme is to fulfill the partial differential algebraic system F(i9t<t), 9s4*, 0) = 0 
at all collocation points and approximate the appearing variables and derivatives at the collocation 
points in terms of (pj, i = 0 ,..., N and j = 0,1,.... We have 

F((9t0,9s0, 0 )^+ 1 %) = 0 , z = 0 ,...,A-l, j = 0 ,l,... 
where ^t^i+i /2 discretized with the following four-point stencils 


0 


7+^ _ ^ 




i-i-1/2 ~ 2 






+‘pr - <pi+i - pd 


i-l-l 
7 _ 

i-l-i 

1-A 


■Pi) 
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Figure A.l. Pattern of the Jacobian for the limit system with m = 9, N = 4 and 
(S) that corresponds to the test case of the 2d Euler-Bernoulli cantilever beam in 
Section EH 


The stencils are based on the convex approximation 0^^'^ = + (1 ~ ^)‘Pi (analogously for 

'^i+ 1 / 2 ) finite differences for the derivative — ipl) /At (analogously for ds’^ij^i/ 2 )- 

This discretization can be in particular viewed as a conservative finite-volume scheme. We obtain 
the unknowns on the new time level for given by solving the system 

= (.Fo(<p^o+Vi+'), = 0 (A.l) 

with = * = 0,...,iV-i. 


Note that in the definition of Ti we treat (9*0,9^0, function in according to 

the discretization stencils stated above. Moreover, the function Q prescribes the boundary conditions 
and (p^ contains the initial conditions. 

In case of the Cosserat rod and Kirchhoff beam models, the function T is nonlinear and (lA.ll) is 
solved with a Newton method. The Jacobi matrix d^pT S K"*(^+i)x"*(^+i) ig a band matrix with 
block structure, presupposing an appropriate sorting of the boundary conditions. The blocks are 
given by 




i-l-1/2'' 


+ Ab -p ^9<,c(0^+^/,) 


As 


where A, B, 9<t>c G '^mxm g^and for the system matrices of the e-dependent model and the limit 
model (e = 0), respectively (cf. (12.8I) - (I2.101) 1. The blocks associated with the limit model are less 
occupied since = Af°^ -l-e^Af^l (analogously for B^, Ce). Moreover, the Jacobian of the discretized 
limit model equals the linear system matrix associated with the discretized model of the first-order 
correction. See Figure lA.ll for an example pattern. 


Remark 9. The presented approach can be also interpreted as spatial semi-discretization with an 
initial value problem as well as temporal semi-discretization with a boundary value problem. In 
any interpretation a differential-algebraic system (DAE) is obtained. Its index varies from 0 to 
3, depending on the underlying system (e > 0 (Cosserat rod), e = 0 (Kirchhoff beam), first-order 
correction 0{e^)), the model variant (M), (T), (S) and the considered semi-discretization, see Ta- 
ble \A.l[ In this context the used discretization schemes can he understood as implicit Runge-Kutta 
methods of stage s = 1.' Gauss method for the DAE in space and in time with X = 0.5 and Radau 
Ila method (Euler method) for the DAE in time with A = 1. Although theoretical convergence results 
are only available up to index 1 for the Gauss method (convergence order p = 2) and up to index 2 
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Table A.l. DAE index for the semi-discretized systems (cf. Remark IH]) 



DAE in time 
(M) (T) (S) 

DAE in space 
(M) (T) (S) 

e > 0 

1 

0 

1 

1 

1 

0 

e = 0 

3 

2 

2 

1 

1 

0 

Oie^) 

3 

2 

2 

1 

1 

0 


for the Radau Ila method (convergence order p = 1) [17] , we point out that in practice the numerical 
scheme can be successfully applied to all cases. 

A.2. Convergence properties. The scheme has the freedom that it can be switched between 
energy-conserving and dissipative by changing the time integration between the Gauss midpoint rule 
(A = 0.5) and the implicit Euler method (A = 1). This change effects the convergence properties, as 
we will show exemplarily for the test case in Section lOI We study at first separately the spatial and 
temporal convergence orders in the context of semi-discretized DAEs (see Remark [9]). As expected 
and in agreement with the analytical results, the midpoint discretization for the spatial DAE yields 
a numerical convergence of order ps = 2 for all variables in both, e-dependent and limit systems. In 
particular, the results are independent of the chosen model variant (M), (T), (S). Figure |A2] (top) 
shows the relative T^(0, l)-error between the reference solution associated with Asref = 10“^ and 
the approximations for As S {10“^, 0.5^ ■ 10“^,..., 0.5^ • 10“^} at T = 2 computed with At = 10“^, 
(S) and A = 1. In contrast, the temporal convergence studies are more sophisticated. For the 
e-dependent system the temporal DAE is of index 0 (T) or 1 (M), (S). In agreement with the 
theory we obtain a numerical convergence of order pt = 2 for A = 0.5 and pt = 1 for A = 1. For 
the limit system the temporal DAE is of index 2 (T), (S) or even 3 (M). In spite of the lack of 
convergence theory we find here pt 1, pt < 1. The algebraic variables behave in general worse 
than the differential ones. Interesting to note is that the results are similar for all model variants. 
Figure [A]2] (bottom) shows the relative T^(0, l)-error between the reference solution associated with 
Atref = 10“^ and the approximations for At G {10“^, 0.5^ • 10“^,..., 0.5’’ • 10“’} at T = 2 computed 
with As = 2 • 10“^, (S) and both A = 0.5 and A = 1. On first glance the error values (in particular 
for the limit system) maybe look large, but this is only due to the used time T. Moreover, we 
point out that the temporal convergence requires a sufficiently small As. Combining temporal 
and spatial convergence and considering At = As —>■ 0, we always obtain the convergence order 
p = 1 for the whole scheme. The only exception is the e-dependent case with moderate e values 
solved with A = 0.5, here we have p = 2. This coincides with the separate investigations as we get 
p = min{pt,Ps}. 
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